Map impacts for all impacted species
Compiled on Mon Aug 17 16:00:43 2020 by ohara
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
source(here('common_fxns.R'))
library(animation)
reload <- FALSE1 Summary
Create a species richness map for threatened marine species.
Combine taxa-level impact maps to create maps of number of threatened species impacted per cell, per year. This is also done by range priority (see the get_spp_priority() function for details)
Plot maps as animated rasters, both by absolute numbers of species impacted in a location and by proportion of species (impact counts / species richness).
This is all done at the level of cumulative impact category, not individual stressors.
2 Methods
2.1 Create species richness map for this assessment
Set up the dataframe of species to include. By taxon, load all species range maps, smash down to taxon-level species richness; then combine these for a total species richness.
spp_incl <- get_incl_spp() %>%
mutate(taxon = tolower(rank))
spp_sens <- spp_incl %>%
filter(!is.na(stressor))
# spp_incl$iucn_sid %>% n_distinct()
# spp_incl %>% filter(!is.na(stressor)) %>% .$iucn_sid %>% n_distinct()
### 1036 if we limit to just sensitive spp, 1271 if we include all
### spp that are mapped/threatened/comp assessedcell_id_rast <- raster(here('_spatial/cell_id_mol.tif'))
ocean_a_rast <- raster(here('_spatial/ocean_area_mol.tif'))
spp_rich_mapfile <- here('_output/rasters/n_spp_map.tif')
if(!file.exists(spp_rich_mapfile) | reload) {
spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
taxa_gps <- spp_incl$taxon %>% unique() %>% sort()
taxa_map_list <- vector('list', length = length(taxa_gps))
for(i in seq_along(taxa_gps)) {
# i <- 2
taxa_gp <- taxa_gps[i]
spp_ids_taxa <- spp_incl %>%
filter(taxon == taxa_gp) %>%
.$iucn_sid %>%
unique()
map_list <- parallel::mclapply(spp_ids_taxa, mc.cores = 12,
FUN = function(id) {
# id <- spp_ids_taxa[1]
spp_f <- file.path(spp_range_dir, sprintf('iucn_sid_%s.csv', id))
spp_df <- read_csv(spp_f, col_types = 'ii') %>%
mutate(iucn_sid = id) %>%
mutate(present = (presence != 5) & !is.na(presence))
})
map_df <- bind_rows(map_list)
map_nspp_df <- map_df %>%
group_by(cell_id) %>%
summarize(n_spp = sum(present))
taxa_map_list[[i]] <- map_to_rast(map_nspp_df, cell_val = 'n_spp')
}
taxa_map_stack <- stack(taxa_map_list)
spp_rich_rast <- raster::calc(taxa_map_stack,
fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(spp_rich_rast, spp_rich_mapfile, overwrite = TRUE)
}
spp_rich_rast <- raster(spp_rich_mapfile)
plot(log(spp_rich_rast),
axes = F, col = hcl.colors(n = 20),
main = 'log(Species richness)')
2.2 Setup for mapping across all taxa
For each impact category (including all):
- gather all taxa-level map files into a raster stack
- collapse the stack using
calc(fun = sum, na.rm = TRUE) - animate over the years
spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
impact_map_dir <- here('_output/rasters/impact_maps')
taxon_impact_map_dir <- file.path(dir_bd_anx, 'taxon_impact_rasts')
all_maps <- list.files(taxon_impact_map_dir,
pattern = '.tif',
full.names = TRUE)
taxa_map_df <- data.frame(map = all_maps) %>%
mutate(year = str_extract(basename(map), '[0-9]{4}') %>% as.integer(),
impact = str_replace_all(basename(map), '.+[0-9]{4}_|.tif', ''))make_ramp <- function(n, palette) {
if(n > 101) {
n_even <- round(n + 50, -2)
breaks <- seq(.5, n_even, length.out = 100)
labels <- seq(0, n_even, by = 100)
} else {
n_even <- n
breaks <- seq(min(0.5, n / 100), n, length.out = 100)
labels <- seq(0, n, by = (n / 10))
}
colors <- hcl.colors(length(breaks), palette = palette)
### For breaks and colors, add a different color for
### a value of exactly zero
return(list('colors' = c('grey40', colors),
'breaks' = c(0, breaks),
'labels' = labels))
}
make_gifs <- function(map_stack, filename, layer_names = NULL) {
if(is.null(layer_names)) layer_names = names(map_stack)
ramp <- make_ramp(max(maxValue(map_stack)), palette = 'viridis')
capture.output({
saveGIF({
for(i in 1:nlayers(map_stack)){
plot(map_stack[[i]],
col = ramp$colors,
breaks = ramp$breaks,
axes = FALSE,
axis.args = list(at = ramp$labels),
main = layer_names[i])
}},
interval = 0.5, movie.name = filename,
ani.width = 700, ani.height = 420)
})
return(invisible(NULL))
}
nspp_rast <- spp_rich_rast
values(nspp_rast)[values(nspp_rast) == 0] <- NA2.3 Create impact map for all included spp, by stressor category
str_cats <- c('fishing', 'climate', 'land-based', 'ocean', 'all')
yrs <- 2003:2013
outfile_stem <- file.path(impact_map_dir, 'impact_%s_%s.tif')
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
impact_map_files <- taxa_map_df %>%
filter(impact == str_cat)
for(y in yrs) {
# y <- 2003
outfile <- sprintf(outfile_stem, str_cat, y)
# message('Processing ', basename(outfile))
if(!file.exists(outfile) | reload) {
imp_yr_files <- impact_map_files %>%
filter(year == y) %>%
.$map
impact_yr_stack <- stack(imp_yr_files)
impact_yr_rast <- raster::calc(impact_yr_stack,
fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(impact_yr_rast, outfile, overwrite = TRUE)
} #else {
# message('File exists: ', outfile, '... skipping!')
# }
}
}outfile <- file.path(impact_map_dir, 'impact_all_2013_2plus.tif')
if(!file.exists(outfile) | reload) {
impact_map_files <- taxa_map_df %>%
filter(impact == 'all_2plus')
impact_2plus_stack <- stack(impact_map_files$map)
impact_2plus_rast <- raster::calc(impact_2plus_stack,
fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(impact_2plus_rast, outfile, overwrite = TRUE)
}2.3.1 maps by percent of threatened species impacted
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
gif_file <- here(sprintf('figs/impact_pct_all_spp_%s.gif', str_cat))
if(!file.exists(gif_file) | reload) {
rast_files <- sprintf(outfile_stem, str_cat, yrs)
map_stack <- stack(rast_files)
x <- map_stack / nspp_rast
if(any(values(x) > 1, na.rm = TRUE)) {
stop('pct raster values greater than 1, wtf')
}
values(x)[values(x) > 1] <- 1
x <- x * 100
make_gifs(x,
filename = gif_file,
layer_names = sprintf('Impacted spp (percent): %s (%s)', str_cat, yrs))
}
# knitr::include_graphics(gif_file)
cat(sprintf('', gif_file))
}
2.3.2 Plot 2013 layers
str_cats <- c('fishing', 'climate', 'land-based', 'ocean', 'all')
outfile_stem <- file.path(impact_map_dir, 'impact_%s_%s.tif')
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
r <- raster(sprintf(outfile_stem, str_cat, 2013))
# r_pri <- raster(sprintf(outfile2_stem, str_cat, 2013))
r_norm <- r / nspp_rast
# r_pri_norm <- r_pri / priority_sum_rast
# plot(r, col = hcl.colors(n = 20),
# axes = FALSE,
# main = paste('count spp impacted:', str_cat, 'stressors'))
plot(r_norm, col = hcl.colors(n = 20),
axes = FALSE,
main = paste('pct spp impacted:', str_cat, 'stressors'))
# plot(r_pri_norm, col = hcl.colors(n = 20),
# axes = FALSE,
# main = paste('priority pct spp impacted:', str_cat, 'stressors'))
}
r <- raster(sprintf(outfile_stem, 'all', '2013_2plus'))
r_norm <- r / nspp_rast
plot(r_norm, col = hcl.colors(n = 20),
axes = FALSE,
main = paste('pct spp impacted: all stressors (2+)'))
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
r_stem <- here('_output/rasters/impact_maps/impact_%s_2013.tif')
imp_r <- raster(sprintf(r_stem, str_cat))
df <- data.frame(impact = values(imp_r),
nspp = values(nspp_rast)) %>%
filter(!is.na(nspp))
df1 <- sample_frac(df, .02)
thingplot <- ggplot(df1, aes(x = nspp, y = impact)) +
geom_point(alpha = .2) +
geom_abline(slope = 1, intercept = 0, color = 'red') +
geom_smooth(method = 'lm', color = 'green', alpha = .5) +
labs(title = sprintf('%s impacts vs spp richness (2 pct sample)', str_cat))
print(thingplot)
}